In [2]:
## input: 1. summarizedExperiement(SE) obj for chromVAR 2. Jaspar matrix
## output: 1. motif x cell (z score) 2. plot: ranked
source("./libs.R")
In [5]:
##------------------------------------------------------------
## inputs
##------------------------------------------------------------
input.chromVar.res.list <- readRDS(file = "../dat/output.jaspar.dev.res.Rdata")
input.chromVar.jaspar.z <- assays(input.chromVar.res.list$dev)$z
input.umap.res <- fread('../dat/1908/Islet_123.MNN_corrected.cluster_labels.filt.txt',header = T)
input.chromVar.jaspar.var <- fread("../dat/1908/output.jaspar.var.res.abcd.csv")
In [6]:
options(repr.plot.width = 2, repr.plot.height = 2)
# filter unkonwn
input.umap.res <- input.umap.res %>% separate(cluster, into = c("cell_type_overall",
"subtype"), remove = F)
head(input.umap.res, 1)
dim(input.umap.res)
ggplot(input.umap.res, aes(UMAP1, UMAP2)) + geom_point(aes(color = cluster), size = 0.1,
alpha = 0.5) + theme_pubr() + theme(legend.position = "none")
In [9]:
require(venn)
options(repr.plot.width = 2, repr.plot.height = 2)
input.umap.res.old <- fread("../dat/Islet_123.MNN_corrected.UMAP.txt")
input.chromVar.jaspar.z <- assays(input.chromVar.res.list$dev)$z
venn(list(old = input.umap.res.old$barcodes, new = input.umap.res$barcodes, motif = colnames(input.chromVar.jaspar.z)[-1]))
table(input.umap.res %>% filter(barcodes %in% colnames(input.chromVar.jaspar.z)[-1]) %>%
pull(cell_type_overall))
table(input.umap.res %>% pull(cell_type_overall))
In [10]:
fun.plot.project.motif <- function(motif, input.chromVar.z = input.chromVar.jaspar.z,
umap.res = input.umap.res, rescale = F, cls, bks = c(-4, 0, 4), ...) {
require(scales)
motif.idx <- grep(motif, rownames(input.chromVar.z))
if (length(motif.idx) == 0) {
message(motif, " is not found!")
return()
}
motif.z <- input.chromVar.z[motif.idx[1], ]
if (rescale) {
sc <- max(abs(quantile(motif.z, probs = c(0.05, 0.95))))
motif.z[motif.z > sc] <- sc
motif.z[motif.z < -sc] <- -sc
}
motif.z <- motif.z %>% as.data.frame() %>% rownames_to_column("barcodes")
colnames(motif.z)[2] <- "zval"
tmp <- umap.res %>% right_join(motif.z)
tmp <- tmp %>% mutate(zval = ifelse(zval > bks[3], bks[3], ifelse(zval < bks[1],
bks[1], zval)))
p.default.cluster.motif <- ggplot(tmp, aes(UMAP1, UMAP2)) + geom_point(aes(colour = zval),
shape = 16, ...) + ggtitle(rownames(input.chromVar.jaspar.z)[motif.idx]) +
scale_color_gradientn(colours = cls, breaks = bks) + theme_pubr() + theme(text = element_blank(),
axis.ticks = element_blank(), legend.position = c(0.2, 0.9), legend.direction = "horizontal",
legend.key.width = unit(2, "mm"), legend.key.height = unit(0.1, "inches"),
legend.text = element_text(size = 10, family = "Arial"), legend.background = element_rect(fill = "transparent",
colour = "transparent"), plot.margin = unit(c(1, 1, -1, -1), "mm"))
p.default.cluster.motif
}
In [12]:
table(input.umap.res%>%pull(cluster))
cell-type specific motif between beta and delta¶
In [14]:
options(repr.plot.width = 4, repr.plot.height = 2)
ps <- lapply(c("PDX1", "FOXA1"), fun.plot.project.motif, bks = c(-2, 0,
2),size=.25,cls=rev(brewer.pal(n = 11, "RdBu")))
names(ps) <- c("PDX1", "FOXA1")
ggarrange(plotlist = ps, ncol = 2)
In [19]:
input.umap.res[is.na(input.umap.res)] <- 0
T Test (one vs other)¶
T Tests (binary)¶
In [15]:
input.chromVar.jaspar.z <- data.table(assays(input.chromVar.res.list$dev)$z,keep.rownames = T)
# aggregate data --------------------------------------------------------------
# melt
input.chromVar.jaspar.z.agg <- melt(input.chromVar.jaspar.z, id = "rn", variable.name = "barcodes",
value.name = "zval")
# add celltype
input.chromVar.jaspar.z.agg <- merge(input.chromVar.jaspar.z.agg, input.umap.res)
head(input.chromVar.jaspar.z.agg,1)
dim(input.chromVar.jaspar.z.agg)
dim(input.chromVar.jaspar.z)
table(input.chromVar.jaspar.z.agg%>%pull(cell_type_overall))
In [184]:
input.chromVar.jaspar.z.agg <- input.chromVar.jaspar.z.agg[complete.cases(input.chromVar.jaspar.z.agg),
] %>% separate(rn, into = c("id", "name"), sep = "_")
range(input.chromVar.jaspar.z.agg$zval)
dim(input.chromVar.jaspar.z.agg)
head(input.chromVar.jaspar.z.agg,1)
In [157]:
table(input.chromVar.jaspar.z.agg %>% filter(name == motif)%>%pull(cell_type_overall))
table(input.chromVar.jaspar.z.agg%>%pull(cell_type_overall))
In [202]:
require(parallel)
celltype.test.all <- list(alpha_vs_gamma = c("alpha", "gamma"), beta_vs_delta = c("beta",
"delta"))
ttest.res <- do.call(rbind, lapply(names(celltype.test.all), function(ntest) {
# cells <- c('beta', 'delta')
celltype.test <- celltype.test.all[[ntest]]
test.motifs <- input.chromVar.jaspar.z.agg %>% filter(cell_type_overall %in%
celltype.test) %>% pull(name) %>% unique
ttest.res <- do.call(rbind, mclapply(test.motifs, function(motif) {
pd <- input.chromVar.jaspar.z.agg %>% filter(name == motif, cell_type_overall %in%
celltype.test)
test.res <- t.test(pd %>% filter(cell_type_overall == celltype.test[1]) %>%
select(zval), pd %>% filter(cell_type_overall == celltype.test[2]) %>%
select(zval))
(data.frame(`motif` = motif, mean_x = test.res$estimate[1], mean_y = test.res$estimate[2],
pval = test.res$p.value/2))
}, mc.cores = 10)) %>% mutate(test = ntest)
}))
head(ttest.res,1)
head(ttest.res%>%filter(test=='beta_vs_delta'),1)
In [333]:
ttest.res <- ttest.res %>% group_by(test) %>% mutate(FDR = p.adjust(pval, "BH"),
padj = p.adjust(pval, "bonferroni")) %>% rowid_to_column("rank") %>% mutate(rank = ifelse(rank >
386, rank - 386, rank)) %>% left_join(input.chromVar.jaspar.z.agg %>% select(name,
id) %>% unique, by = c(motif = "name"))
In [387]:
fwrite(ttest.res, "~/Dropbox (UCSD_Epigenomics)/workReports/2019-08-27_islet/chromVar_ttest_res.csv")
In [16]:
ttest.res<- fread("~/Dropbox (UCSD_Epigenomics)/workReports/2019-08-27_islet/chromVar_ttest_res.csv")
dim(ttest.res)
Motif db annotate¶
install.packages('data.tree')
In [48]:
head(df,1)
In [95]:
df <- ttest.res.2 %>% filter(FDR < 0.01) %>% select(family.id, contains("name")) %>%
unique() %>% drop_na %>% mutate(name = "all", family.id = paste0("0.", family.id)) %>%
unite(path,name,name.y, name.x, family.name, sep = ";")
df %>% head(2)
dt.tmp <- FromDataFrameTable(df, pathName = "path", pathDelimiter = ";")
dt.tmp
In [140]:
df <- ttest.res.2 %>% filter(FDR < 0.01) %>% select(contains("id"), contains("name")) %>%
unique() %>% drop_na %>% mutate(name = "all", family.id = paste0("0.", family.id))
df %>% head(1)
In [163]:
data(MisLinks)
data(MisNodes)
str(MisLinks)
str(MisNodes)
In [168]:
MisNodes.2 <- rbind(df %>% select(starts_with("family")) %>% rename(id = "family.id",
name = "family.name"), df %>% select(class.id, name.x) %>% rename(id = "class.id",
name = "name.x"), df %>% select(superclass.id, name.y) %>% rename(id = "superclass.id",
name = "name.y"), data.frame(id = "-1", name = "all")) %>% rowid_to_column("rowid") %>%
mutate(id = factor(id))
MisNodes.2 %>% str
MisNodes.dic <- MisNodes.2$rowid-1
names(MisNodes.dic) <- MisNodes.2$name
In [169]:
MisLinks.2 <- ToDataFrameNetwork(dt.tmp)%>%mutate(value=1)
MisLinks.2$from = MisNodes.dic[MisLinks.2$from]
MisLinks.2$to = MisNodes.dic[MisLinks.2$to]
MisLinks.2%>%str
In [170]:
# Load data
data(MisLinks)
data(MisNodes)
# Plot
forceNetwork(Links = MisLinks.2, Nodes = MisNodes.2,
Source = "from", Target = "to",
Value = "value", NodeID = "name",
Group = "id", opacity = 0.8)
In [143]:
head(ToDataFrameNetwork(dt.tmp),1)
In [114]:
diagonalNetwork( ToListExplicit(dt.tmp, unname = T),fontSize = 10)
In [96]:
simpleNetwork(ToDataFrameNetwork(dt.tmp), fontSize = 6,)
In [124]:
require(igraph)
net <- graph.data.frame(links, nodes, directed=T)
In [122]:
install.packages(c("network",'sna'))
In [121]:
require(GGally)
require(network)
#require(sna)
net = rgraph(10, mode = "graph", tprob = 0.5)
net = network(net, directed = FALSE)
In [ ]:
In [88]:
acmeNetwork
In [87]:
require(networkD3)
acmeNetwork <- ToDataFrameNetwork(acme, "name")
simpleNetwork(acmeNetwork[-3], fontSize = 12)
In [136]:
options(repr.plot.width = 12, repr.plot.height = 12)
dt.tmp.ig <- as.igraph(dt.tmp, directed = TRUE, direction = "climb")
#net <- simplify(dt.tmp.ig, remove.multiple = F, remove.loops = T)
plot(dt.tmp.ig,vertex.size=6,edge.arrow.size=.25)
In [116]:
require(igraph)
plot(as.igraph(dt.tmp, directed = TRUE, direction = "climb"),cex=.5)
In [82]:
SetGraphStyle(dt.tmp,rankdir='TB')
plot(dt.tmp)
In [42]:
require(dendextend)
dt.tmp %>% as.dendrogram%>%set('branches_k_color',k=1)%>%plot
In [22]:
dim(ttest.res)
ttest.res.2 <- (ttest.res %>% left_join(tfclass.db.dic$merged %>% select(jaspar.id, family.id,
family.name), by = c(id = "jaspar.id"))%>%unique())
In [55]:
head(tfclass.db$class,1)
In [61]:
ttest.res.2<- ttest.res.2 %>% mutate(class.id = sub(".[[:digit:]]$", "", family.id)) %>%
left_join(tfclass.db$class%>%select(-about), by = c("class.id" = "id"))%>%
mutate(superclass.id = sub(".[[:digit:]]$", "", class.id)) %>%
left_join(tfclass.db$superclass%>%select(-about), by = c("superclass.id" = "id"))
In [514]:
tmp.dup <- ttest.res.2 %>% group_by(test) %>% filter(duplicated(motif))
#ttest.res.2 %>% filter(motif %in% tmp.dup$motif)
In [400]:
sum(!ttest.res$id %in% tfclass.db.dic$jasparTOtfclass$jaspar.id)
sum(!ttest.res$id %in% tfclass.db.dic$merged$jaspar.id)
remain.motif <- ttest.res[!ttest.res$id %in% tfclass.db.dic$merged$jaspar.id,c('motif','id')]%>%unique
sum(remain.motif$motif %in% tfclass.db.dic$merged$genus.name)
remain.motif[!(remain.motif$motif %in% tfclass.db.dic$merged$genus.name),]
In [17]:
tfclass.db.dic <- readRDS("~/github/atacMotif/db/dic_jaspar_tfclass.rds")
str(tfclass.db.dic)
In [57]:
tfclass.db <- readRDS('~/github/atacMotif/db/tfclass.rds')
str(tfclass.db)
In [360]:
rm(input.pseudotime)
plots - ranked bars¶
In [2]:
#require(devtools)
#install_github("jokergoo/ComplexHeatmap")
require(ComplexHeatmap)
In [331]:
ttest.res <- ttest.res %>% mutate(padj.1 = ifelse(enriched %in% c("gamma", "delta"),
-FDR, FDR))%>%arrange(enriched,padj.1)%>%select(-padj.1)
In [301]:
ttest.res <- ttest.res %>% mutate(diff = mean_x - mean_y) %>% arrange(test, diff)
In [17]:
require(circlize)
In [296]:
options(repr.plot.width = 2.5, repr.plot.height = 4)
ttest.res %>% head(1)
pd.rank <- ttest.res %>% filter(test == "alpha_vs_gamma") %>% ungroup
# pd.rank <- ttest.res %>% filter(test == 'beta_vs_delta')
require(ComplexHeatmap)
pvalue_col_fun = colorRamp2(c(0, 2, 3, 5), c("white", "grey100", "grey50", "grey20"))
ha = rowAnnotation(enriched = anno_simple(pd.rank %>% pull(enriched), col = c(alpha = "darkred",
gamma = "magenta4")), padj = anno_simple(-log10(pd.rank %>% pull(padj)), col = pvalue_col_fun))
ht = Heatmap(pd.rank %>% select(starts_with("mean")) %>% rename(alpha = mean_x, gamma = mean_y) %>%
as.matrix(), right_annotation = ha, cluster_rows = F, name = "mean", cluster_columns = F)
lgd_pvalue = Legend(title = "padj", col = pvalue_col_fun, at = c(0, 2, 3, 5), labels = c("1",
"0.01", "0.001", "0.00001"))
lgd_enriched = Legend(title = "enriched", legend_gp = gpar(fill = c("darkred", "magenta4")),
labels = c("alpha", "gamma"))
draw(ht, annotation_legend_list = list(lgd_pvalue, lgd_enriched))
In [310]:
options(repr.plot.width = 2.5, repr.plot.height = 4)
ttest.res %>% head(1)
# pd.rank <- ttest.res %>% filter(test == 'alpha_vs_gamma') %>% ungroup
pd.rank <- ttest.res %>% filter(test == "beta_vs_delta") %>% ungroup %>% filter(padj <
0.01)
require(ComplexHeatmap)
pvalue_col_fun = colorRamp2(c(0, 2, 3, 50), c("white", "grey100", "grey50", "grey20"))
cls <- coolwarm(7)
main_col_fun = colorRamp2(c(-2, -0.5, 0, 0.5, 2), cls[c(1:2, 4, 6:7)])
ha = rowAnnotation(enriched = anno_simple(pd.rank %>% pull(enriched), col = c(beta = "green",
delta = "gold")), padj = anno_simple(-log10(pd.rank %>% pull(padj)), col = pvalue_col_fun),
diff = anno_simple(pd.rank$diff, col = colorRamp2(c(-1.3, 0, 1.1), c("green",
"white", "gold"))))
ht = Heatmap(col = main_col_fun, pd.rank %>% select(starts_with("mean")) %>% rename(beta = mean_x,
delta = mean_y) %>% as.matrix(), right_annotation = ha, cluster_rows = F, name = "mean",
cluster_columns = F)
lgd_pvalue = Legend(title = "padj", col = pvalue_col_fun, at = c(0, 2, 3, 50), labels = c("1",
"0.01", "0.001", "1e-50"))
lgd_enriched = Legend(title = "enriched", legend_gp = gpar(fill = c("green", "gold")),
labels = c("beta", "delta"))
draw(ht, annotation_legend_list = list(lgd_pvalue, lgd_enriched))
In [353]:
pd.rank%>%head(20)%>%mutate(rank=motif.1)%>%select(rank,motif)
In [350]:
options(repr.plot.width = 3, repr.plot.height = 4)
require(pals)
pd.rank <- ttest.res %>% filter(test == "beta_vs_delta") %>% ungroup
head(pd.rank, 1)
n_m =nrow(pd.rank %>% filter(enriched == 'delta'))
n_m =0
pd.rank <- pd.rank %>% mutate(motif.1 = as.numeric(factor(motif, levels = pd.rank$motif))-n_m,
lpadj = ifelse(enriched == "beta", log10(FDR), -log10(FDR))) %>% filter(FDR <
0.01) %>% mutate(enriched = factor(enriched, levels = c("delta", "beta")))
ggplot(pd.rank) + geom_bar(aes(motif.1, lpadj, fill = abs(lpadj)), stat = "identity") +
theme_pubr(base_size = 10) + theme() + coord_flip(expand = F) + scale_fill_gradientn(colours = pal.compress(coolwarm)[7:13]) +
scale_x_continuous(breaks = seq(0-n_m, 386-n_m, 20)) + facet_grid(enriched ~ ., scales = "free_y",
space = "free_y")+geom_hline(yintercept = 0,linetype=2)
# geom_vline(xintercept = %>%
# pull(rank))) +
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.b_d.rank.pdf",
width = 3, height = 4, units = "in")
In [349]:
pd.rank%>%head(10)
In [347]:
options(repr.plot.width = 3, repr.plot.height = 4)
require(pals)
pd.rank <- ttest.res %>% filter(test == "alpha_vs_gamma") %>% ungroup
head(pd.rank, 1)
n_m = nrow(pd.rank %>% filter(enriched == "gamma"))
n_m = 0
pd.rank <- pd.rank %>% mutate(motif.1 = as.numeric(factor(motif, levels = pd.rank$motif)) -
n_m, lpadj = ifelse(enriched == "alpha", log10(padj), -log10(padj))) %>% filter(padj <
0.01) %>% mutate(enriched = factor(enriched, levels = c("gamma", "alpha")))
ggplot(pd.rank) + geom_bar(aes(motif.1, lpadj, fill = abs(lpadj)), stat = "identity") +
theme_pubr(base_size = 10) + theme() + coord_flip(expand = F) + scale_fill_gradientn(colours = pal.compress(coolwarm)[7:13],
breaks = c(0, 2, 5, 20, 40, 60, 80)) + scale_x_continuous(breaks = seq(0 -
n_m, 386 - n_m, 20)) + facet_grid(enriched ~ ., scales = "free_y", space = "free_y") +
geom_hline(yintercept = 0, linetype = 2)
# geom_vline(xintercept = %>% pull(rank))) +
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.a_g.rank.pdf",
width = 3, height = 4, units = "in")
plots¶
In [223]:
cols.subtypes=c(alpha1='darkred',alpha2='lightpink',
beta1='darkgreen',beta2='lightgreen',
gamma1='orange',gamma2='gold')
In [237]:
cols.celltype["beta"] = rgb(t((col2rgb(cols.celltype["beta_1"]) + col2rgb(cols.celltype["beta_2"]))/2),
maxColorValue = 255)
cols.celltype["alpha"] = rgb(t((col2rgb(cols.celltype["alpha_1"]) + col2rgb(cols.celltype["alpha_2"]))/2),
maxColorValue = 255)
cols.celltype["delta"] = rgb(t((col2rgb(cols.celltype["delta_1"]) + col2rgb(cols.celltype["delta_2"]))/2),
maxColorValue = 255)
In [185]:
options(repr.plot.width = 3, repr.plot.height = 2)
m <- "NFIL3"
celltype.test <- c("alpha", "gamma")
plotViolin <- function(motif = m, test_celltypes = celltype.test) {
pd <- input.chromVar.jaspar.z.agg %>% filter(name == motif, cell_type_overall %in%
test_celltypes)
ggviolin(pd, x = "cell_type_overall", remove = T, width = 0.5, y = "zval", size = 0.5,
shape = 16, draw_quantiles = c(0.25, 0.5, 0.75), fill = "cell_type_overall",
ylab = "Motif enrichment") + theme_pubr() + ggtitle(motif) + geom_hline(yintercept = 0,
linetype = 2, size = 0.25) + theme(legend.position = "none", axis.text.x = element_text(angle = 45,
hjust = 1), axis.title.x = element_blank())
}
options(repr.plot.width = 1.5, repr.plot.height = 2)
plotViolin("SCRT1", c("beta", "delta"))
In [346]:
plotViolin("SCRT1", c("beta", "delta")) + scale_fill_manual(values = c(beta = cols.celltype$beta,
delta = cols.celltype$delta))
plotViolin("TCF7L2", c("beta", "delta")) + scale_fill_manual(values = c(beta = cols.celltype$beta,
delta = cols.celltype$delta))
In [253]:
ttest.res%>%filter(motif %in% c("TCF7L2"))
In [145]:
options(repr.plot.width = 4, repr.plot.height = 4)
p1<- plotMotif(m="MAFG", "alpha_vs_gamma",size=1,bks=c(-3,0,3))
p2<-plotMotif(m="MEIS1", "alpha_vs_gamma",size=1,bks=c(-3,0,3))
#plotMotif(m="ELF5", "alpha_vs_gamma",size=1,bks=c(-2,0,2))
plotMotif(m="EHF", "alpha_vs_gamma",size=1,bks=c(-2,0,2))
In [356]:
options(repr.plot.width = 2, repr.plot.height = 2)
plotMotif(m="FOXP2", "alpha_vs_gamma",size=.5, cls = rev(vals),bks=c(-2,0,2))
In [355]:
options(repr.plot.width = 2, repr.plot.height = 2)
plotMotif(m="FOXP2", "beta_vs_delta",size=.5, cls = rev(vals),bks=c(-2,0,2))
In [345]:
plotMotif(m="TCF7L2", "alpha_vs_gamma",size=.5, cls = rev(vals),bks=c(-2,0,2))
In [159]:
options(repr.plot.width = 2, repr.plot.height = 2)
plotMotif <- function(m = "SCRT1", cmp = "alpha_vs_gamma",...) {
p <- fun.plot.project.motif(motif = m,...)
p %+% (p$data %>% filter(cell_type_overall %in% unlist(strsplit(cmp, split = "_vs_"))))
}
plotMotif(m="MAFK", "alpha_vs_gamma",size=.5)
In [339]:
pd.rank%>%filter(grepl('TCF',motif))
In [340]:
options(repr.plot.width = 2, repr.plot.height = 2)
p <- plotMotif("TCF3", "beta_vs_delta", size = 0.75, bks = c(-3, 0, 3), cls = rev(vals))
p
In [329]:
vals <- brewer_pal(palette = "RdBu")(7)
vals[3:5] <- vals[4]
options(repr.plot.width = 2, repr.plot.height = 2)
p <- plotMotif("EHF", "beta_vs_delta", size = 0.75, bks = c(-3, 0, 3), cls = rev(vals))
p
p <- plotMotif("EHF", "beta_vs_delta", size = 0.75, bks = c(-2, 0, 2), cls = rev(vals))
p
p <- plotMotif("EHF", "beta_vs_delta", size = 0.5, bks = c(-2, 0, 2), cls = rev(vals))
p
In [255]:
options(repr.plot.width = 2, repr.plot.height = 2)
p <- plotMotif("EHF", "beta_vs_delta", size = 0.5, bks = c(-3, 0, 3))
xxmin = -1.3
xxmax = 1
yymin = -2.5
yymax = 1.5
p + geom_rect(xmin = xxmin, xmax = xxmax, ymin = yymin, ymax = yymax, fill = NA,
color = "black")
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.ehf.b_d.umap.pdf",
width = 2, height = 2, units = "in")
p$data <- p$data %>% filter(cell_type_overall == "delta")
p + coord_cartesian(ylim = c(yymin, yymax), xlim = c(xxmin, xxmax)) + theme(legend.position = "none")
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.ehf.b_d.umap.insert.pdf",
width = 1, height = 1, units = "in")
In [247]:
options(repr.plot.width = 1.5, repr.plot.height = 2)
plotViolin("EHF", c("beta", "delta")) + scale_fill_manual(values = c(beta = cols.celltype$beta,
delta = cols.celltype$delta))
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.ehf.b_d.umap.insert.violin.pdf",
width = 1.5, height = 2, units = "in")
In [250]:
plotMotif("MAFF", "beta_vs_delta", size = 0.5, bks = c(-2, 0, 2))
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.maff.b_d.umap.pdf",
width = 2, height = 2, units = "in",useDingbats=F)
In [249]:
options(repr.plot.width = 1.5, repr.plot.height = 2)
plotViolin("MAFF", c("beta", "delta")) + scale_fill_manual(values = c(beta = cols.celltype$beta,
delta = cols.celltype$delta))
ggsave("~/Dropbox (UCSD_Epigenomics)/workReports/2019-09_islet_rev/Fig.1.maff.b_d.umap.insert.violin.pdf",
width = 1.5, height = 2, units = "in")
In [208]:
#options(repr.plot.width = 6, repr.plot.height = 6)
p<- plotMotif("EHF", "beta_vs_delta",size=.5,bks=c(-2,0,2),alpha=.75)
p
plotViolin("EHF", c("beta", "delta"))
plotMotif("MAFF", "beta_vs_delta",size=.5,bks=c(-3,0,3),alpha=.75)
plotViolin("MAFF", c("beta", "delta"))
In [341]:
options(repr.plot.width = 2, repr.plot.height = 2)
plotMotif("EGR", "alpha_vs_gamma")
In [339]:
options(repr.plot.width = 2, repr.plot.height = 2)
plotMotif("EHF", "beta_vs_delta")
In [334]:
options(repr.plot.width = 5, repr.plot.height = 4, max_z = 3)
plotSNmotif <- function(testt = "alpha_vs_gamma", fdr_th = 0.001) {
top_motifs <- ttest.res %>% filter(padj < fdr_th & test == testt) %>% pull(motif)
cell_types <- unlist(strsplit(testt, "_vs_"))
pd <- input.chromVar.jaspar.z.agg %>% filter(name %in% top_motifs & cell_type_overall %in%
cell_types) %>% select(barcodes, name, zval, cell_type_overall) %>% mutate(motif = as.numeric(factor(name,
levels = top_motifs)))
dim(pd)
pd$zval[pd$zval > max_z] <- max_z
pd$zval[pd$zval < -max_z] <- -max_z
tmp <- pd %>% select(barcodes, cell_type_overall) %>% unique() %>% arrange(cell_type_overall)
# dim(tmp %>% mutate(barcodes))
pd <- pd %>% mutate(barcodes = as.numeric(factor(barcodes, levels = tmp$barcodes)))
print(seps <- table(tmp %>% unique() %>% pull(cell_type_overall)))
print(seps.motif <- table(ttest.res %>% filter(padj < fdr_th & test == testt) %>%
pull(enriched)))
ncolor = 5
ggplot(pd, aes(x = barcodes, y = motif)) + geom_raster(aes(fill = zval)) + geom_hline(yintercept = as.numeric(cumsum(seps.motif))[1],
linetype = 2, size = 0.5) + geom_vline(xintercept = as.numeric(cumsum(seps))[1],
linetype = 2, size = 0.5) + coord_cartesian(expand = F) + theme_bw() + theme(panel.grid = element_blank()) +
scale_fill_gradientn(colors = colorRampPalette(c("blue", "white", "red"))(ncolor),
breaks = seq(-3, 3, length.out = ncolor))
}
In [337]:
plotSNmotif(testt = 'alpha_vs_gamma',fdr_th = .01)
In [377]:
p <- plotSNmotif(testt = "alpha_vs_gamma", fdr_th = 0.01)
In [400]:
quantileNorm <- function(dat.vect,lowq = 0.05) {
dat <- as.numeric(dat.vect)
lowq = 0.05
(qs <- as.numeric(quantile(dat, c(lowq, 1 - lowq))))
dat <- (dat - qs[1])/(qs[2] - qs[1])
dat[dat > 1] <- 1
dat[dat < 0] <- 0
dat
}
In [424]:
require(pals)
In [441]:
p <- plotSNmotif(testt = "beta_vs_delta", fdr_th = 0.001)
p$data <- p$data %>% group_by(motif) %>% mutate(zval = scale(zval, center = T))
p$data$zval[p$data$zval > 2] <- 2
p$data$zval[p$data$zval < -2] <- -2
print(p + scale_fill_gradientn(colours = colorRampPalette((viridis(30)), bias = .75)(30)))
In [429]:
p <- plotSNmotif(testt = "alpha_vs_gamma", fdr_th = 0.01)
p$data <- p$data %>% group_by(motif) %>% mutate(zval = scale(zval, center = T))
p$data$zval[p$data$zval > 2] <- 2
p$data$zval[p$data$zval < -2] <- -2
print(p + scale_fill_gradientn(colours = colorRampPalette(parula(30),bias=.75)(30)))
In [406]:
options(repr.plot.width = 5, repr.plot.height = 4)
p <- plotSNmotif(testt = "alpha_vs_gamma", fdr_th = 0.01)
p$data <- p$data %>% group_by(motif) %>% mutate(zval = quantileNorm(dat.vect = zval))
print(p + scale_fill_gradientn(colours = colorRampPalette(cols.hm.avg.tf(30), bias = 0.5)(30),
breaks = seq(0, 1, length.out = 3)))
In [348]:
options(repr.plot.width = 5, repr.plot.height = 4, max_z = 3)
p<- plotSNmotif(testt = 'alpha_vs_gamma',fdr_th = .01)
print(p)
In [335]:
plotSNmotif(testt = 'beta_vs_delta')
In [ ]: